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1 Introduction 


Black holes have always been at the centre of theoretical interests in quan¬ 
tum gravity and string theory. Through the non-perturbative formulation 
of string theory and quantum gravity via gauge/string correspondence jTj, 
significant progresses have been made in our understanding of black hole 
microstates, especially its entropy counting. However, we still need to un¬ 
derstand better various aspects of quantum gravitational system, and one of 
the important aspects is understanding their real-time evolution. 

One model, which is more tractable to analyse the real-time dynamics 
compared with many other theoretical models in the gauge/string corre¬ 
spondence, is a BFSS (Banks-Fischler-Shenker-Susskind) matrix model [2]. 
This model is a 0 + 1 dimensional large N matrix quantum mechanics, and 
conjectured to be dual to either ten-dimensional string or eleven di¬ 

mensional M theory depending on the various parameter limit. Since BFSS 
matrix model is dual to the higher dimensional spacetime, it would be very 
interesting to ask how the higher dimensional black hole formation can be 
understood from the dual matrix model. In this set-up, we consider the 
various head-on collisions of two bunches of DO-branes and their real-time 
evolution, which is expected to become a single bound state of DO-branes at 
late time. This is a matrix model thermalization and dual to the formation 
of one large black hole in either ten- or eleven dimension. Our main interest 
in this paper is, by solving the matrix model classical equation of motion, 
to study quantitatively how such a bound state is formed in the real-time 
evolution of the large N matrix model. 

In order to study the complicated time-evolution of large N matrix model, 
we make several simplifications; First, we consider only the classical limit 
in this paper, namely we solve the classical equation of motion numerically 
and find its time-dependent solution. At first thought, this simplification 
might sound very boring limit. However, there are several reasons why clas¬ 
sical limit of late time evolution is interesting and capture aspects of non- 
perturbative physics. Even though classical limit is justified in the weak 
coupling limit, since what we conduct is not the perturbative analysis near 
the trivial vacuum but rather the analysis seeking the time-dependent soli- 
ton configuration, this can capture non-perturbative physics. Remember that 
thermalization is a very complicated but also universal phenomenon which 
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one cannot see in various solvable models, nor in perturbation analysis from 
the trivial vacuum [5J. It has a property that all memories of the initial data 
are lost in late time, and this late time thermalization may occur only in the 
large N limit. Therefore even though we simply solve the classical equations 
of motion, we are looking some non-perturbative physics through the late 
time evolution. 

Another simplification is, that we neglect effects of fermions for simplicity. 
Since supersymmetry is expected to be more crucial at low energy but not in 
high energy deconfined phase, we expect that, for late time thermalization, 
the analysis with/without supersymmetry does not change much. In our 
analysis, we have conduct N up to N = 24, and regard it as large enough to 
perform an extrapolation to the large N limit. 

Before we close the introduction, we comment on several related references. 
Thermal equilibrium (i.e., static) properties of the BFSS matrix models for 
black holes are very well studied. Thermodynamic simulations of BFSS ma¬ 
trix model was initiated by Kabat et al PE] by the mean-field method [5], 
and recently by Monte Carlo method mmm- These results confirmed that 
the duality is valid not just at supergravity level but also at stringy level; 
at finite-coupling ra and finite-TV m- For non-equilibrium properties, a 
probe limit thermalization was studied in m lUS Ifij- For so-called BMN 
matrix modeQ Berenstein et al studied the thermalization numerically in 
the series of papers [18.1, T9l I20J in detail. BFSS matrix model with initial 
conditions different from ours have also been studied in Finally 

the black hole formulation at the correspondence point for the BFSS matrix 
model was analytically studied recently in [2111221 In this paper we consider 
the similar setting to these works, though the parameter range is different, 
which is possible since we solve the model numerically. 


2 The model and simulation method 

BFSS matrix model [2] is the maximally supersymmetric matrix quantum 
mechanics, which is the dimensional reduction of 4d J\f = 4 supersymmetric 

2 BMN matrix model m is a model which has mass for the adjoint scalars. This model 
is supposed to be dual to the 11-dimensional pp-waves. 
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Yang-Mills theory to one dimension (time). In this paper, we concentrate on 
the classical dynamics, and set fermion to zero for simplicity. Henceforth, we 
consider only the bosonic part of the Lagrangian of BFSS matrix model and 
it is given by 

L = I Y(AW,) 2 + \ Y [V„, V,] 2 ) . (1) 

\ m=1 z m ^ n ) 

Here Xm (M = 1, 2, • • • ,9) are N x N Hermitian matrices and (DiXm) 
is the covariant derivative given by ( D t XM ) = where A t 

is the U(N) gauge field. This matrix model is a gravity-decoupled theory 
on the N DO-branes [23]. Large N bound state of DO-branes with finite 
temperature T with large’t Hooft coupling is dual to a black hole (black 0- 
brane) with Hawking temperature T [3j . Such a bound state is described by 
highly non-commutative matrices. On the other hand, if matrices are close 
to block-diagonal, it describes multi-DO-branes and each of which is highly 
quantum)^] black hole. 

In the A t = 0 gauge, the classical dynamics is described by the equation of 
motion 

^XL~YyX N ,[X M ,X N ]] = 0 (2) 

N 

with the Gauss’s law constraint 

= 0. (3) 

In the rest of this paper we solve this classical equation of motion. From 
the action ([TJ it is clear that ’t Hooft coupling A = gy M N appears only as 
an overall factor. Therefore the classical limit, i.e., h —» 0 limit, is equivalent 
to the effective coupling constant A e // = X/U 3 goes to zero limit in the ’t 
Hooft scaling limit, where N —y oo and U is the typical VEV scale in the 
matrix model associated with the temperature, or quantum fluctuations of 
the matrix degrees of freedom. Therefore, our analysis of solving classical 
equation of motion is justified in hight temperature/energy limit^J 

2 In order to have classical black hole picture we need both large’t Hooft coupling and 
large N. For single DO-brane, JV = 1 therefore highly quantum description. 

3 From the scaling symmetry t t/a, Xm —> cxXm of classical equation of motion Q, 


E 


Xmi 


dX 


M 


dt 
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2.1 Discretization 


In order to study the time evolution, we discretize the equation of motion 
0 while preserving the constraint @ exactly. For that purpose, we write 
the EOM as 


T , U\ dx M(t) „ dV M (t) 

v u(t)= dt , F M (t)= dt 

(4) 

and 


F M (t) = [X M (t),X N (t)}}. 

j 

(5) 

The constraint becomes 


J2l x M(t),V M (t)}= 0. 

(6) 


M 


Under the infinitesimal time evolution t —> t + St, X M and Vu change as 
X M (t + St) = X M (t) + V M (t) • St + F M (t ) • + O ((St) 3 ) , (7) 

V M (t + St) = V M (t) + (F M (t ) + F M (t + St)) • | + 0 ((St) 3 ) , (8) 

as one can easily check by using the standard Taylor expansion. We terminate 
this expansion at the order of (St) 2 , 

X M (t + St) = X M (t) + V M (t) • St + F M (t) • (9) 

V M (t + St) = V M (t) + (F M (t) + F M (t + St)) ■ |. (10) 

Then, a simple but tedious calculation shows that the constraint 0 is sat¬ 
isfied at t + St if it holds at t, without any additional higher-order terms: 

(t),V M (t)]=0 =+ (t + St ), V M (t + St)] — 0. (ii) 

M M 

one might wonder if we can re-scale the total energy. However this is a symmetry only at 
the classical level and full quantum theory does not possess such a symmetry. 
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Therefore, we use i§ 
Q at t = 0, 


and (10), and choose the initial configuration to satisfy 


^[X M (0),y M (0)]=0. (12) 

M 

2.2 Initial condition — collision of two bunches of D0- 
branes — 

We consider a collision of two bunches of DO-branes, which is analogous to 
two black zero-branes, or two black holes (the left of Fig. [Tj) . The initial 
condition i,^ 



V 2 = --- = Vg = 0, (13) 

and 

= (x j ;r = +^,«) (i4) 

for n — 2, 3 , • • • , 9, i = 1,2,--- , N/2 and j = N/2+1, ■ ■ ■ ,N, where a^ij and 
are Gaussian random numbers with the weights e - “ 2 / 2 and e _&2//2 . Here 
a is the source for off-diagonal elements and therefore non-commutativity, 
which is crucial ingredients for the bound state formation. If a = 0 , then 
commutator square potential vanishes and two bunches just pass by. For 
nonzero a, the quantitative details of the initial condition can depend on 
the choice of random numbers, especially when N is small. Note that we 
consider only even values of N. This initial condition satisfies the Gauss- 
law constraint Y1m=i[Xm,Vm] = 0 . We consider the’t Hooft large -N limit, 
where A = Qy M N is fixed, L and V are fixed, and then, the energy scales as 
N 2 . When two bunches are separated, the off-diagonal element a has larger 

4 It would be better to introduce the off-diagonal elements also in the diagonal blocks. 
For simplicity, we do not do that for the moment. (If we introduce them, then the width 
should be taken larger than that for the off-diagonal blocks, a.) 
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mass as a and so quantum mechanically path-integrated out perturbatively 
in X/a 3 and suppressed. This describes the quantum fluctuation of open 
strings stretched between two bunches. Because the potential energy at t = 
0 is approximately ^Tt[Xm, X^] 2 ~ N ■ N 2 L 2 a 2 , the natural scaling is 
a ~ l/{yNL). We will see that the energy actually scales as E ~ N 2 then. 


3 Formation of a large quantum black hole 

3.1 Formation of a single bound state 

Let us first ask ‘when a single bound state of eigenvalues is formed’. This is 
the process where two bunches of DO-branes, or equivalently two quantum 
black holes merge to form one large quantum black hole (Fig. [Tj). 



Figure 1: If a single bound state of eigenvalues is formed, it is analogous to a 
formation of one large black hole, and this is the dual to the thermalization 
process of the whole matrix system. 

The easiest way to see this process is to plot TrXf^, especially TvX 2 , as 
a function of t. In Fig. i TvXf/N for N = 8, for L = 5.0, V = 0 and 
y/Ncr = 0.12, with several different values of time step dt, is shown. Here 
we set the ’t Hooft scale A = 1. At late time, errors associated with discrete 
time steps become non-negligible. We can see that behaviors at t < 70 can 
reliably be studied with dt = 0.00010. TrXf/N bounces a few times and 
then stays small, which suggests the formation of a single bound state. Note 
that, if instead two bunches oscillated as shown in Fig. I TrX'l/N would 
oscillate without decaying. 

At late time, we expect that the system goes to the the “typical” configu¬ 
rations and that such typical states are rotationally symmetric after taking 
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the average over time and/or different initial configurations. In Fig. [4j the 
average (TrXj^/N) (M = 1, 2, • • • , 9), by using 0 < t < 10, 10 < t < 20, • • •, 
60 < t < 70, for 100 different initial configurations, are plotted. We can see 
that the average values converge to the same value at late time, suggesting 
the restoration of the rotational symmetry. 


3.1.1 Insensitivity to initial conditions for thermalization 

In order to add further evidence for thermalization, we studied two com¬ 
pletely different initial conditions with the same value of the energy. In 
order to tune the total energy, we set Vm = 0 and rescale X M to aX m- Then 
the energy scales as E —>■ a 4 E. By using this, we can tune the energy to 
any value. Here we consider (i) N = 4, L = 10.0, V = 0, VNa = 0.12, 
and (ii) N = 4, all matrix components are generated Gaussian weight. Then 
we rescale Xm so that E/N 2 = 0.5. Note that (i) is much more anisotropic 
initial condition than (ii). We employed dt = 0.00001 and dt = 0.000001 
in this calculation, which give consistent results at t < 90. From Fig. [5| we 
can see that X^A 4 -=i TrJ 5 T^ f /9A7 with two different initial conditions behave 
the same manner in late time when the energies are the same, and it is hard 
to tell the initial condition unless we explicitly solve the equation of motion 
to go back to the past. These suggest typical thermliazation, i.e., whatever 
initial conditions it starts with, a system ends up with similar typical states. 


3.2 Large -N limit and Thermalization 

In order to study the statistical nature of the BFSS matrix model, it is not 
essential to take the step size dt very small, as long as the energy is conserved 
and we take enough late time evolution. In the following, we take dt = 0.0005. 


3.2.1 Large- N limit with fixed L, V and \/No 

Let us first consider a large- N limit with fixed L. V and y/Ncr. As a concrete 
example, we consider L = 5.0, V = 0 and \/Na = 0.12. Similar results were 
obtained for other values too. 
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In Fig. [6| we plot Tr Xf/N for N = 8,12,16. We can see qualitatively 
similar behaviours, and the fluctuation at late-time becomes smaller as N 
becomes larger. When a bound state is formed, the Virial theorem relates 
the kinetic energy K and the potential energy V as (K) = 2(V) = | E, where 
( • ) stands for the time average and E = K + V is the total energy, which is 
conserved. Therefore, (K — | E)/N 2 is a good measure for the fluctuation. 
As we can see from Fig.[7j this quantity fluctuates around zero and suppressed 
more as N becomes larger. 

In order to see the statistical property at late-time, we take the time average 
at 50 < t < 100 and then take an average over random initial configurations 
(50 samples for N = 4,8, 20 samples for N = 12, 15 samples for N = 16 
and 10 samples for N = 24). The error bar is estimated purely statistically. 
The results are plotted in Fig. [8| especially the bottom figure of Fig. [8] shows 
how N dependence appears in (| K — | E\). Note that since | K — ‘/E\ is an 
absolute value, its fluctuation is always added up. The fact that 

lim (|K - ?B|) = 0 (15) 

iV—>-oo 3 

suggests that the fluctuation is suppressed as N becomes larger. In finite 
N , (l-h'— ^E\/N 2 ^ is proportional to 1/N, while the corrections to E/N 2 
and (Y^m=i TrW^/OlV) are proportional to 1/W 2 . E/N 2 is consistent with 
ci + C 2 /N 2 behaviour with some constant c\ and C 2 , where C 2 could be zero. 


3.2.2 Large- N limit with fixed E/N 2 


Let us consider another, perhaps more natural large -N limit, in which E/N 2 
is fixed exactljj^} We first generate initial configurations with L = 5.0, V = 0 
and \f~Na = 0.12, and then set E/N 2 to 1.5, by rescaling scalars as we did in 
the end of §|3.1[ We collected 50 samples for N = 4,8, 20 samples for N = 12, 


15 samples for N = 16 and 10 samples for N = 24. We use 50 < t < 100 
again. In principle, since the energy is fixed, the time average for all samples 
with different initial configurations should give the same average value within 
errors for a sufficiently long time interval. (Still, we terminate the simulation 
at t = 100 so that the conservation of the energy is not violated due to 


5 In the large N limit, which corresponds to the thermodynamic limit, this is the same 
as fixing the temperature of the system. 
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the discretization error.) Therefore, statistical fluctuations of the average 
(Em=i^m/^) become much smaller compared to Fig.JsJ The results 
are plotted in Fig. fel We can see that the correction to (Y1m =i TrXf^/QA") 
also starts with 1/N 2 . 


3.2.3 Large-N limit and thermalization 

As we have seen, the system shows typical thermalization process and the 
fluctuation of macroscopic quantities disappears at large N. Therefore, the 
macroscopic quantities converge to the same value at late time, irrespectively 
of the initial condition, as long as the energy is the same. The information 
of the initial condition is hidden in the 1/N correction. This is exactly what 
we expect. Since the large N limit corresponds to the thermodynamic limit 
there is no distinguish between micro canonical ensemble (this is what we 
did) and canonical ensemble where temperature is fixed. The difference only 
appears in 1/A" expansion. Note that an exact thermalization occurs only in 
the large N limit. 


3.3 Thermalization time 


Let us investigate the thermalization time. There are several different ways to 
define the thermalization time. Although it involves the dual interpretation, 
perhaps one of the most natural ones is the time scale for the “black hole” 
formation. In Fig. 10, we plot the history of (TvXf/N) for N = 8,12,16, L = 
5.0, V = 1.0 and yN a = 0.16. We used 10 different sets of random numbers 
and took average. Data at N = 12 and N = 16 are indistinguishable within 
errors, suggesting that the time evolution of (TrXf/N) has a well-defined 
large -N limit. Remember that large N limit is a thermodynamic limit, where 
canonical ensemble and microcanonical ensemble are distinguishable only by 
looking at the 1/N suppressed order. Therefore, we call the system get 
thermalized at the late-time if the system becomes indistinguishable “typical 
states”, and the deviations for the macroscopic quantities from the typical 
states are only of order 1/N. Then, to read-off the thermalization time, i.e., 
the time-scale system get thermalized, we have to read-off the decay of the 
deviation at late-time in great detail to the order of 1/N. However this is 
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hard generically since the error bar is too large and we have initial-condition- 
dependence for the late-time average. 


To avoid this issue, we define the thermalization time by another way. 
This is defined by the time scale for a small perturbation added on top of 
the bound state to become invisible. Note that if N is not sufficiently large, 
then the thermalization time defined in this way would depend on the detail 
of the perturbation. Let us consider two-point function along time direction 
which does not require the perturbation, and from that, we will read off the 
thermalization time. Let us consider a gauge-invariant operator 

—Tv(XM(s)W s ^ + tX M (s + t)Wg t 

s+t)? (16) 

where W s . s+t is the Wilson fine, 


W. 




= Pe*W~ dt 'Xt') 


(17) 


In our setup we took the A = 0 gauge, so 

^Tr (X M (s)W s , s+t X M (s + t)Wl +t ) = ^Tr (X M (s)X M (s + t)). (18) 

Let us consider 

m = Y, ( 4tV(V«(s)V„( S + «))), (19) 

M ' ' 


where the time average ( • ) is taken at late-time. We have already seen that 
/(0) is of order one. When t becomes large, f(t) decays because Xm(s + 
t ) “forgets” the information of Xm{s). Therefore, it should be possible to 
determine the thermalization time from the time scale for the decay of /(£). 
The numerical result is shown in Fig. 11 It turned out that f(t) can be 
fitted well as 


f(t) ~ ce r * 2 cos (cot ), (20) 

where the overall constant c = (yfTrXjvf(f) 2 ), the decay width T and 
the frequency c o are of order N°. The width T naturally gives the time scale 
for the thermalization, 

(Thermalization time) 1/VT~N°. (21) 
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We can also consider connected two-point functions of gauge invariant op¬ 
erators. We consider one of the operators studied by Berenstein and col¬ 
laborators, TV (XM(t)Xnr(t)) (M ^ N ). Although we did not find a simple 
fitting function for this, the decay is consistent with e~ r 1 (Fig. 12). Then 
the thermalization time defined by this operator is 


(Thermalization time) ~ 1/T 7 ~ N°. (22) 

Fig. [12] suggests that the thermalization time, defined from the exponential 
decay, converges to finite value in the large N limit. This is quite interesting 
since it is consistent with finiteness of quasinormal modes in the black hole 
background, although this convergence is more subtle at late time due to 
finite N effect. As is discussed in [24], although the correlation function keep 
decaying to zero at any finite order in the 1/W-expansion, the decay should 
disappear once full finite-IV corrections are taken into account. The fat tail 
at long time scale in Fig. [12] would demonstrate this property. 


3.4 Under what initial conditions, does thermalization 
occur? 

Let us consider at what values of V and a the thermalization, or equivalently 
the formation of the bound state, can be realized. For simplicity we fix L to 
be 5.0. Firstly let us note that the off-diagonal elements are the source for 
the non-linearity, which is responsible for the formation of the bound state. 
If we set the off-diagonal elements to be zero (a = 0) at t = 0, then the off- 
diagonal elements are not generated by the classical equations of motion and 
no interaction appears at all. In such case, two bunches just pass through 
each other. 

Even if a is nonzero, if it is too small, the non-linearity cannot grow large 
enough before the collision. Therefore, it is clear that a must be sufficiently 
large for the bound state formation. However, as we will see shortly, it turns 
out that a should not be too small but also must not be too large in order 
to form the bound state. This indicates that if the a accelerates the two 
bunches before the first collision too much, they tend to pass by without 
forming the bound state. 
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In order to illustrate these, we show the evolution of TrX? [N for N = 16, 
for L = 5.0, V = 1.0 and various values of y/Na in Fig. [l3j For y/Na = 
0.16,0.12 and 0.08, the late-time behaviours at t > 20 look qualitatively 
the same. A fact that a bounce at y/N a = 0.12 is larger than those at 
y/No = 0.08 and y/Na = 0.16 suggests very complicated dynamics. At 
y/N a = 0.04, on the other hand, two bunches just pass through with each 
other. It might be possible that the bunches comes back later after very long 
time and form a bound state; in fact sometimes we observe that TV X^/N 
becomes as large as 0(100) and then come back and form a single bound 
state. As we can see from the early-time behavior in Fig. [13, 


as a increases 


the bunches are accelerated more before the collision (i.e. Tr Xf/N decreases 
more rapidly). So whether a single bound state corresponding to thermalised 
state is formed or not depends on a very subtle competition of the acceleration 
before the collision and deceleration after the collision, both of which are 
enhanced by a larger a. 



bouncing too much. We can see that a single bound state is not formed when 
y/N a is too large, which means the acceleration before the collision beats the 
deceleration after the collision. We cannot get any conclusive pattern from 
these, probably the result depends on the choice of the random numbers. We 
left this issue as open questions. 


4 Summary and discussions 


In this paper, we considered the head-on collisions of two bunches of D0- 
branes and study its real-time evolution in the BFSS matrix model. Espe¬ 
cially our interests in this study are how time evolution differs and under 
what conditions the big bound states corresponding thermalized state are 
formed and what is the its time-scale. 

We have seen that all of TrA | f /N where M = I, • • • ,9 converges at late 
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time to the same values and rotational symmetry is restored irrespective to 
the anisotropic initial conditions in §3.1[ Especially the insensitivity to the 


initial conditions strongly suggest the thermalization and formation of a large 


bound state dual to quantum black holes in late time. In §3.2[ for various 
scaling limit, we have checked that at the thermalized stage, Virial theorem, 
which is expected to hold at the equilibrium state, is in fact satisfied and 
deviation from that occurs only at the subleasing order in 1/N expansion. 


In 1T3, we estimated thermalization time from the time-direction two point 
functions and its exponential decay at late time, and shows that the timescale 
is N° order, which seems to converge in the large N limit. This result is 
very interesting since it is consistent with quasi-normal mode in the dual 
black hole. Finally in §3.4[ we studied the initial condition dependence for 
the thermalization. Clearly for small <x, the off-diagonal sources for non¬ 
commutativity, two DO-brane bunches pass by while we notice that for too 
large cr, two bunches also pass by in such cases. We however have no clear 
interpretation for this results. 

Our analysis is conducted under several assumptions. The big assumptions 
we take in our analysis is that we can completely neglect the quantum effects 
and fermion effects. Quantum effects are important as the effective coupling 
becomes larger, and fermions are also very important since they change the 
interactions between two bunches of DO-branes. It is quite possible that once 


we take into account these effects, then the time-scale we estimated at §3.3 
can be modified significantly. Obviously more detail study is necessary to 
get conclusive results. 


To better understand under what initial conditions does the system ther¬ 
malized, let us discuss what is the typical states in the BFSS matrix model. 
If the final states is made up by sets of almost commuting matrices, then 
the entropy of the gas is made up by N DO-branes, which has only O(N) 
entropy. Off-diagonal elements are highly suppressed in such case and does 
not contribute to the entropy. On the other hand, a single large black hole, 
which is a single bound state of the DO-branes, has the entropy of order 
N 2 , because all of the off-diagonal elements are excited. We can also con¬ 
sider multi-black hole states. Clearly a single black hole is entropically and 
therefore free-energy-wise dominant, since all of the off-diagonal elements are 
equally excited and contributes to the entropy. 


13 







However there are subtle issues. In the BFSS matrix model, because moduli 
space is not bounded, when the two bunches pass through each other, it is 
possible, at least classically, that they just pass by and run to infinity. This 
is all due to the fact that all of the adjoint scalars have zero masfj^J However 
by taking into account quantum fluctuations, this is unlikely, because in 
the real system involving quantum fluctuation, there are always small but 
nonzero quantum fluctuation and when the non-linearity from the interaction 
is large enough at late time, two bunches are likely to be pulled each other 
and merge to a single bound state. If the entropy of the black hole is large 
enough (i.e. N is large enough), then it is unlikely to see a large deviation 
from the black hole within a finite simulation time. 

Note that since our calculation is classical, it is crucial to take small ef¬ 
fective coupling limit i.e., A e ff —»■ 0 to justify our results. Even though 
effecting coupling is small, the equation of motion can accumulate large non¬ 
linearity at late time, which can capture non-perturbative soliton formation 
physics. A large N effect (which plays the same role as thermodynamic limit) 
and the late time universality effects (which capture enough non-linearity of 
the interaction and scramble) are crucial ingredients for the thermalization. 
There is an argument that the classical time evolution of the BFSS matrix 
model is stochastic [25]. Such natures are probably crucial for the system 
get thermalized and show the universal behaviour at late time. 
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Figure 2: TrXf/N for N = 8, for L = 5.0, V = 0 and y/Na = 0.12; 
dt = 0.00025,0.00010 and 0.00001. The horizontal axis is time t. We can 
see that dt = 0.00025 gives correct answer up to t ~ 60. Beyond there, the 
error becomes large quickly, even for dt = 0.00010. The bottom panel is a 
zoom-up of the top one. 
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Figure 3: A ‘bound state’ would not necessarily be a single black hole; for 
example, it would be possible that two bunches oscillates at the late time. 
However such a configuration cannot be regarded as thermalized state. 
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Figure 4: ( TrX^/N ) (M = 1,2, ••• ,9), N = 8, for L = 5.0, V = 0 and 
dt = 0.00010, t/Nct = 0.12 (top) and \/N<j = 0.14 (bottom). The values 
along the horizontal axis are the time range used for taking the time average. 
The error bar is purely statistical. We can see the restoration of the rotational 
symmetry at late time. 
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Figure 5: TrX^/OfV, N = 4, E/N 2 = 0.5, dt = 0.000001. Blue 

and red lines represent the initial conditions (i) and (ii) in the main text, 
respectively. The bottom panel is the zoom up of the top one. The late time 
behaviors are indistinguishable. 
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Figure 6: Y^m=i TrX^/9iV for N = 8,12 and 16, for L = 5.0, V = 0 
and y/N<r = 0.12. The bottom panel is a zoom-up of the top one. Qualita¬ 
tively similar behaviours can be seen and the fluctuation at late-time becomes 
smaller as N becomes larger. 
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Figure 7: —jj%— for iV = 8,12 and 16, for L = 5.0, V = 0 and \fNa = 0.12. 
Here K and E are the kinetic energy and total energy. When the bound 

_ 2 jjj 

state is formed, the long-time average of N % must be zero due to the virial 
theorem. The fluctuation around zero becomes smaller as N becomes larger. 
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Figure 8: ^"-dependence of (X)m=i TrX^/91V) (top), E/N 2 (middle) and 
(\K — | E\/N 2 ) (bottom). Collision parameters are L = 5.0, V = 0 and 
VNa = 0.12. 24 
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Figure 9: iV-dependence of (Em=i (top) and (\K — ^E\/N 2 } 

(bottom), where the energy at each run is fixed to E/N 2 = 1.5. 

(Eki Tr ^/9iV) behaves as const. + const./iV 2 , while (| K - \E\/N 2 ) is 
consistent with const. /N behaviour. 
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Figure 10: One-point function (TrXf/N) for N = 8,12,16, L = 5.0, V = 1.0 
and y/Na = 0.16. Average over 10 samples. Results with different values of 
N are not distinguishable within statistical error. 
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Figure 11: f(t) = J2 m <^Tt(X m (s)X m (s + 1))) for N = 6,8,16, and 24, 
L = 5.0, V — 0 and VNc t = 0.16. The fitting line is for N = 24, with the 
ansatz fit) ~ c e~ rt2 cos (cut) and fit parameters c = 0.182, T = 0.0311 and 
uj = 1.48. Better convergence to N = oo can be seen at earlier time. 
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Figure 12: Logarithm of ^ ((TrX II (s)X u (s))(TrX IJj (s + t)X„(s + t))) as a 
function of t, for N = 8,16, and 24, L = 5.0, V = 0 and VlSfcr = 0.16. 
The fitting line is for N = 24, with 0.062 • e -a245 L Better convergence to 
N = oo can be seen at earlier time. 
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Figure 13: TrX^ for N = 16, for L = 5.0, V = 1.0 and various values of 
VN<T. 
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Figure 14: N = IQ, L = 5.0. The green “x” symbols show (' V,y/Na ) in 
which two bunches passed through with each other and Tr Xf/N became 
at least twice bigger than the initial value. A single-bunch state would be 
formed after long time, but not immediately. The red “+” symbols show 
(V, y/Ncr) where a single bound state is formed without bouncing too much. 
We used the same set of random numbers ay and by for all V and \fNa. 
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